import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import rasterio as rio
from rasterio.plot import show
from rasterio.plot import show_hist
import os
import fiona
!pwd
/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis
!ls
Approx_SG.kml debug.log Displaying images.ipynb Extracting_SG_from_Landsat8_images.ipynb Generated DataFrames HDB data Landsat_Image_processing_final_experiments.ipynb Landsat8_image_processing_initial_experimentation.ipynb Population_data_World_Bank.ipynb Population_density_and_Housing_density_analysis.ipynb Predicting the future of SG.ipynb SG_cropped_Landsat8_images SG_Landsat8_images SG_shapefiles
os.chdir('SG_Landsat8_images/')
!pwd
/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis/SG_Landsat8_images
!ls
LC08_L1TP_125059_20130627_20170503_01_T1.tar LC08_L1TP_125059_20140427_20170423_01_T1.tar LC08_L1TP_125059_20141020_20170418_01_T1.tar LC08_L1TP_125059_20150601_20170408_01_T1.tar LC08_L1TP_125059_20160923_20170320_01_T1.tar LC08_L1TP_125059_20170419_20170501_01_T1.tar LC08_L1TP_125059_20180929_20181009_01_T1.tar LC08_L1TP_125059_20190916_20190925_01_T1.tar LC08_L1TP_125059_20200310_20200314_01_T1.tar
os.chdir('LC08_L1TP_125059_20200310_20200314_01_T1.tar/')
!ls
Bands_2020 = {}
for i in range(1,12):
var = 'imgB'+str(i)
img_file = 'LC08_L1TP_125059_20200310_20200314_01_T1_B' + str(i) + '.TIF'
Bands_2020[var] = rio.open(img_file, 'r')
#print(Bands_2020)
os.chdir('..')
LC08_L1TP_125059_20200310_20200314_01_T1_ANG.txt LC08_L1TP_125059_20200310_20200314_01_T1_B1.TIF LC08_L1TP_125059_20200310_20200314_01_T1_B10.TIF LC08_L1TP_125059_20200310_20200314_01_T1_B11.TIF LC08_L1TP_125059_20200310_20200314_01_T1_B2.TIF LC08_L1TP_125059_20200310_20200314_01_T1_B3.TIF LC08_L1TP_125059_20200310_20200314_01_T1_B4.TIF LC08_L1TP_125059_20200310_20200314_01_T1_B5.TIF LC08_L1TP_125059_20200310_20200314_01_T1_B6.TIF LC08_L1TP_125059_20200310_20200314_01_T1_B7.TIF LC08_L1TP_125059_20200310_20200314_01_T1_B8.TIF LC08_L1TP_125059_20200310_20200314_01_T1_B9.TIF LC08_L1TP_125059_20200310_20200314_01_T1_BQA.TIF LC08_L1TP_125059_20200310_20200314_01_T1_MTL.txt
os.chdir('LC08_L1TP_125059_20190916_20190925_01_T1.tar')
!ls
Bands_2019 = {}
for i in range(1,12):
var = 'imgB'+str(i)
img_file = 'LC08_L1TP_125059_20190916_20190925_01_T1_B' + str(i) + '.TIF'
Bands_2019[var] = rio.open(img_file, 'r')
#print(Bands_2019)
os.chdir('..')
LC08_L1TP_125059_20190916_20190925_01_T1_ANG.txt LC08_L1TP_125059_20190916_20190925_01_T1_B1.TIF LC08_L1TP_125059_20190916_20190925_01_T1_B10.TIF LC08_L1TP_125059_20190916_20190925_01_T1_B11.TIF LC08_L1TP_125059_20190916_20190925_01_T1_B2.TIF LC08_L1TP_125059_20190916_20190925_01_T1_B3.TIF LC08_L1TP_125059_20190916_20190925_01_T1_B4.TIF LC08_L1TP_125059_20190916_20190925_01_T1_B5.TIF LC08_L1TP_125059_20190916_20190925_01_T1_B6.TIF LC08_L1TP_125059_20190916_20190925_01_T1_B7.TIF LC08_L1TP_125059_20190916_20190925_01_T1_B8.TIF LC08_L1TP_125059_20190916_20190925_01_T1_B9.TIF LC08_L1TP_125059_20190916_20190925_01_T1_BQA.TIF LC08_L1TP_125059_20190916_20190925_01_T1_MTL.txt
os.chdir('LC08_L1TP_125059_20180929_20181009_01_T1.tar')
!ls
Bands_2018 = {}
for i in range(1,12):
var = 'imgB'+str(i)
img_file = 'LC08_L1TP_125059_20180929_20181009_01_T1_B' + str(i) + '.TIF'
Bands_2018[var] = rio.open(img_file, 'r')
#print(Bands_2018)
os.chdir('..')
LC08_L1TP_125059_20180929_20181009_01_T1_ANG.txt LC08_L1TP_125059_20180929_20181009_01_T1_B1.TIF LC08_L1TP_125059_20180929_20181009_01_T1_B10.TIF LC08_L1TP_125059_20180929_20181009_01_T1_B11.TIF LC08_L1TP_125059_20180929_20181009_01_T1_B2.TIF LC08_L1TP_125059_20180929_20181009_01_T1_B3.TIF LC08_L1TP_125059_20180929_20181009_01_T1_B4.TIF LC08_L1TP_125059_20180929_20181009_01_T1_B5.TIF LC08_L1TP_125059_20180929_20181009_01_T1_B6.TIF LC08_L1TP_125059_20180929_20181009_01_T1_B7.TIF LC08_L1TP_125059_20180929_20181009_01_T1_B8.TIF LC08_L1TP_125059_20180929_20181009_01_T1_B9.TIF LC08_L1TP_125059_20180929_20181009_01_T1_BQA.TIF LC08_L1TP_125059_20180929_20181009_01_T1_MTL.txt
os.chdir('LC08_L1TP_125059_20170419_20170501_01_T1.tar')
!ls
Bands_2017 = {}
for i in range(1,12):
var = 'imgB'+str(i)
img_file = 'LC08_L1TP_125059_20170419_20170501_01_T1_B' + str(i) + '.TIF'
Bands_2017[var] = rio.open(img_file, 'r')
#print(Bands_2017)
os.chdir('..')
LC08_L1TP_125059_20170419_20170501_01_T1_ANG.txt LC08_L1TP_125059_20170419_20170501_01_T1_B1.TIF LC08_L1TP_125059_20170419_20170501_01_T1_B10.TIF LC08_L1TP_125059_20170419_20170501_01_T1_B11.TIF LC08_L1TP_125059_20170419_20170501_01_T1_B2.TIF LC08_L1TP_125059_20170419_20170501_01_T1_B3.TIF LC08_L1TP_125059_20170419_20170501_01_T1_B4.TIF LC08_L1TP_125059_20170419_20170501_01_T1_B5.TIF LC08_L1TP_125059_20170419_20170501_01_T1_B6.TIF LC08_L1TP_125059_20170419_20170501_01_T1_B7.TIF LC08_L1TP_125059_20170419_20170501_01_T1_B8.TIF LC08_L1TP_125059_20170419_20170501_01_T1_B9.TIF LC08_L1TP_125059_20170419_20170501_01_T1_BQA.TIF LC08_L1TP_125059_20170419_20170501_01_T1_MTL.txt
os.chdir('LC08_L1TP_125059_20160923_20170320_01_T1.tar')
!ls
Bands_2016 = {}
for i in range(1,12):
var = 'imgB'+str(i)
img_file = 'LC08_L1TP_125059_20160923_20170320_01_T1_B' + str(i) + '.TIF'
Bands_2016[var] = rio.open(img_file, 'r')
#print(Bands_2016)
os.chdir('..')
LC08_L1TP_125059_20160923_20170320_01_T1_ANG.txt LC08_L1TP_125059_20160923_20170320_01_T1_B1.TIF LC08_L1TP_125059_20160923_20170320_01_T1_B10.TIF LC08_L1TP_125059_20160923_20170320_01_T1_B11.TIF LC08_L1TP_125059_20160923_20170320_01_T1_B2.TIF LC08_L1TP_125059_20160923_20170320_01_T1_B3.TIF LC08_L1TP_125059_20160923_20170320_01_T1_B4.TIF LC08_L1TP_125059_20160923_20170320_01_T1_B5.TIF LC08_L1TP_125059_20160923_20170320_01_T1_B6.TIF LC08_L1TP_125059_20160923_20170320_01_T1_B7.TIF LC08_L1TP_125059_20160923_20170320_01_T1_B8.TIF LC08_L1TP_125059_20160923_20170320_01_T1_B9.TIF LC08_L1TP_125059_20160923_20170320_01_T1_BQA.TIF LC08_L1TP_125059_20160923_20170320_01_T1_MTL.txt
!pwd
/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis/SG_Landsat8_images
def show_bands(dictionary):
Bands = dictionary
fig, ax = plt.subplots(3,3, figsize = (15,15))
num = 0
for i in range(3):
for j in range(3):
num += 1
key = 'imgB' + str(num)
ax[i,j].imshow(Bands[key].read(1))
ax[i,j].set_xlabel(key)
def show_histogram(dictionary):
Bands = dictionary
fig, ax = plt.subplots(3,3, figsize = (15,15))
num = 0
for i in range(3):
for j in range(3):
num += 1
key = 'imgB' + str(num)
show_hist(Bands[key], bins = 50, ax = ax[i,j])
#ax[i,j].set_xlabel(key)
!pwd
/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis/SG_Landsat8_images
os.chdir('..')
!pwd
/c/Users/Bharath/anaconda3/envs/SG_Land_Reclamation_analysis
show_bands(Bands_2020)
show_histogram(Bands_2020)
show_bands(Bands_2019)
show_histogram(Bands_2019)
show_bands(Bands_2018)
show_histogram(Bands_2018)
show_bands(Bands_2017)
show_bands(Bands_2016)
fig, ax = plt.subplots(1,2, figsize = (15,15))
ax[0].imshow(Bands_2016['imgB5'].read(1))
ax[1].imshow(Bands_2016['imgB1'].read(1))
<matplotlib.image.AxesImage at 0x2291075ed00>
plt.figure(figsize = (20,10))
plt.imshow((Bands_2016['imgB5'].read(1) - Bands_2016['imgB1'].read(1))[3700:4700, 1500:3700])
<matplotlib.image.AxesImage at 0x2290ef16400>
plt.figure(figsize = (20,10))
plt.imshow((Bands_2016['imgB5'].read(1) - Bands_2016['imgB4'].read(1))[3700:4700, 1500:3700])
<matplotlib.image.AxesImage at 0x2290f12d5e0>
plt.figure(figsize = (20,10))
plt.imshow((Bands_2017['imgB5'].read(1) - Bands_2017['imgB4'].read(1))[3700:4700, 1500:3700])
<matplotlib.image.AxesImage at 0x2293e86c640>
plt.figure(figsize = (20,10))
plt.imshow((Bands_2018['imgB5'].read(1) - Bands_2018['imgB4'].read(1))[3700:4700, 1500:3700])
<matplotlib.image.AxesImage at 0x2291725fa60>
plt.figure(figsize = (20,10))
plt.imshow((Bands_2019['imgB5'].read(1) - Bands_2019['imgB4'].read(1))[3700:4700, 1500:3700])
<matplotlib.image.AxesImage at 0x2293e7d64c0>
(Bands_2019['imgB5'].read(1)).shape
(7741, 7591)
(Bands_2019['imgB9'].read(1)).shape
(7741, 7591)
plt.figure(figsize = (15,15))
plt.imshow(Bands_2019['imgB5'].read(1))
<matplotlib.image.AxesImage at 0x2290ef448e0>
plt.figure(figsize = (15,15))
plt.imshow(Bands_2019['imgB9'].read(1))
<matplotlib.image.AxesImage at 0x22906fdb820>
This is a major challenge for analyzing Landsat images for locations like Singapore that are very cloudy almost throughout the year.
def remove_clouds(dictionary):
Bands = dictionary
img_without_cloud = Bands['imgB5'].read(1) - Bands['imgB9'].read(1)
fig, ax = plt.subplots(1,3, figsize = (15,15))
ax[0].imshow(Bands['imgB5'].read(1))
ax[0].set_xlabel('imgB5')
ax[1].imshow(Bands['imgB9'].read(1))
ax[1].set_xlabel('imgB9')
ax[2].imshow(img_without_cloud)
ax[2].set_xlabel('img_without_cloud')
return img_without_cloud
wc_2020 = remove_clouds(Bands_2020)
plt.imshow(wc_2020)
<matplotlib.image.AxesImage at 0x2293e608430>
plt.figure(figsize=(20,10))
plt.imshow(wc_2020[3700:4700, 1500:3700])
<matplotlib.image.AxesImage at 0x2290ef61490>
wc_2019 = remove_clouds(Bands_2019)
plt.imshow(wc_2019)
<matplotlib.image.AxesImage at 0x229170b0a00>
plt.figure(figsize=(20,10))
plt.imshow(wc_2019[3700:4700, 1500:3700])
<matplotlib.image.AxesImage at 0x22916ff20d0>
wc_2018 = remove_clouds(Bands_2018)
plt.figure(figsize=(20,10))
plt.imshow(wc_2018[3700:4700, 1500:3700])
<matplotlib.image.AxesImage at 0x22906fbce80>
wc_2017 = remove_clouds(Bands_2017)
plt.figure(figsize=(20,10))
plt.imshow(wc_2017[3700:4700, 1500:3700])
<matplotlib.image.AxesImage at 0x2290ef45a00>
wc_2016 = remove_clouds(Bands_2016)
plt.figure(figsize=(20,10))
plt.imshow(wc_2016[3700:4700, 1500:3700])
<matplotlib.image.AxesImage at 0x2290f30c190>
def NormalizedDiff(band1, band2, band9):
b1_wo_clouds = band1.read(1)
b2_wo_clouds = band2.read(1)
#b1_wo_clouds = band1.read(1) - band9.read(1)
#b2_wo_clouds = band2.read(1) - band9.read(1)
b1 = b1_wo_clouds.astype('float')
b2 = b2_wo_clouds.astype('float')
ND = np.where(b1+b2==0., -1.0, (b1 - b2)/(b1 + b2))
return ND
NDWI_2020 = NormalizedDiff(Bands_2020['imgB3'],Bands_2020['imgB6'], Bands_2020['imgB9'] )
<ipython-input-57-ba7b939c5a5e>:12: RuntimeWarning: invalid value encountered in true_divide ND = np.where(b1+b2==0., -1.0, (b1 - b2)/(b1 + b2))
plt.imshow(NDWI_2020)
<matplotlib.image.AxesImage at 0x229170d9ee0>
plt.figure(figsize=(20,10))
plt.imshow(NDWI_2020[3700:4800, 1500:3700])
<matplotlib.image.AxesImage at 0x22916e7b070>
NDWI_2019 = NormalizedDiff(Bands_2019['imgB3'],Bands_2019['imgB6'], Bands_2019['imgB9'])
<ipython-input-57-ba7b939c5a5e>:12: RuntimeWarning: invalid value encountered in true_divide ND = np.where(b1+b2==0., -1.0, (b1 - b2)/(b1 + b2))
plt.imshow(NDWI_2019)
<matplotlib.image.AxesImage at 0x2290f0955e0>
plt.figure(figsize=(20,10))
plt.imshow(NDWI_2019[3700:4700, 1500:3700])
<matplotlib.image.AxesImage at 0x229106ce160>